require(estimatr)

#### data preparation ####
online.data <- read.csv("online_data.csv")

placebo.data <- read.csv("placebo_data.csv")
nrow(placebo.data)  # number of respondents

# dichotomize education
online.data$college <- (online.data$education > 7) * 1
placebo.data$college <- (placebo.data$education > 7) * 1

#### sensitivity analysis ####
DQ.est.online <- mean(online.data$DQ, na.rm = TRUE)
DQ.se.online <- sqrt(DQ.est.online * (1 - DQ.est.online) / 
                       sum(! is.na(online.data$DQ)))

list.result.online <- lm_robust(list ~ treatment, data = online.data, 
                                fixed_effects = ~ block)
list.est.online <- as.numeric(list.result.online$coefficients)
list.se.online <- as.numeric(list.result.online$std.error)

list.result.placebo <- lm_robust(list ~ treatment, data = placebo.data, 
                                 fixed_effects = ~ block)
list.est.placebo <- as.numeric(list.result.placebo$coefficients)
list.se.placebo <- as.numeric(list.result.placebo$std.error)
round(list.est.placebo, 3)
round(list.se.placebo, 3)

# object labels are from Equation (A.3)
lambda <- seq(0.24, 0.99, 0.01)  # percentage of non-strategic misreporters
pi.P <- c(0, 0.05, 0.1)  # true approval rate of Paul Nueva
pi.M <- array(NA, c(76, 3, 6))  # true approval rate of Duterte
for (i in 1:3) {
  corrected.list.est <- as.numeric(list.est.online - list.est.placebo) / 
    (1 - lambda)
  corrected.list.se <- 
    as.numeric(sqrt((list.se.online ^ 2 + list.se.placebo ^ 2) / 
                      (1 - lambda) ^ 2))
  corrected.SDB.se <- 
    as.numeric(sqrt(DQ.se.online ^ 2 + 
                      (list.se.online ^ 2 + list.se.placebo ^ 2) / 
                      (1 - lambda) ^ 2))
  # corrected list estimate
  pi.M[, i, 1] <- pi.P[i] + corrected.list.est
  pi.M[, i, 2] <- pi.P[i] + corrected.list.est + 
    qnorm(0.025) * corrected.list.se
  pi.M[, i, 3] <- pi.P[i] + corrected.list.est + 
    qnorm(0.975) * corrected.list.se
  # corrected SDB estimate
  pi.M[, i, 4] <- DQ.est.online - (pi.P[i] + corrected.list.est)
  pi.M[, i, 5] <- DQ.est.online - (pi.P[i] + corrected.list.est) + 
    qnorm(0.025) * corrected.SDB.se
  pi.M[, i, 6] <- DQ.est.online - (pi.P[i] + corrected.list.est) + 
    qnorm(0.975) * corrected.SDB.se
}

# proportion of non-strategic misreporters when the point estimate of 
# Duterte's approval rate reach its DQ-based estimate
1 - sum(pi.M[, 1, 1] >= DQ.est.online) * 0.01
# proportion of non-strategic misreporters when the upper limit of
# the 95% CI reach the DQ-based estimate
1 - sum(pi.M[, 1, 3] >= DQ.est.online) * 0.01

# SDB estimate and the 95% CI when varying the proportion of
# non-strategic misreporters from 24% to 50%
round(pi.M[24 - 23, 1, 4:6], 3)
round(pi.M[30 - 23, 1, 4:6], 3)
round(pi.M[40 - 23, 1, 4:6], 3)
round(pi.M[50 - 23, 1, 4:6], 3)

#### Figure 5 ####
cairo_pdf("Figure_5.pdf", width = 6, height = 4, pointsize = 8)
layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(4, 4, 1, 1), oma = c(0, 0, 2, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.71, 0.01), seq(0.71, 0.24, -0.01)), 
        c(pi.M[pi.M[, 1, 1] < 1, 1, 2], rev(pi.M[pi.M[, 1, 1] < 1, 1, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M[, 1, 1] < 1], pi.M[pi.M[, 1, 1] < 1, 1, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.7, 0.01), seq(0.7, 0.24, -0.01)), 
        c(pi.M[pi.M[, 2, 1] < 1, 2, 2], rev(pi.M[pi.M[, 2, 1] < 1, 2, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M[, 2, 1] < 1], pi.M[pi.M[, 2, 1] < 1, 2, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.68, 0.01), seq(0.68, 0.24, -0.01)), 
        c(pi.M[pi.M[, 3, 1] < 1, 3, 2], rev(pi.M[pi.M[, 3, 1] < 1, 3, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M[, 3, 1] < 1], pi.M[pi.M[, 3, 1] < 1, 3, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.71, 0.01), seq(0.71, 0.24, -0.01)), 
        c(pi.M[pi.M[, 1, 1] < 1, 1, 5], rev(pi.M[pi.M[, 1, 1] < 1, 1, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M[, 1, 1] < 1], pi.M[pi.M[, 1, 1] < 1, 1, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.7, 0.01), seq(0.7, 0.24, -0.01)), 
        c(pi.M[pi.M[, 2, 1] < 1, 2, 5], rev(pi.M[pi.M[, 2, 1] < 1, 2, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M[, 2, 1] < 1], pi.M[pi.M[, 2, 1] < 1, 2, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.68, 0.01), seq(0.68, 0.24, -0.01)), 
        c(pi.M[pi.M[, 3, 1] < 1, 3, 5], rev(pi.M[pi.M[, 3, 1] < 1, 3, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M[, 3, 1] < 1], pi.M[pi.M[, 3, 1] < 1, 3, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Approval rate of \"Paul Nueva\" = 0", at = 0.19, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.05", at = 0.19 + 1 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.1", at = 0.19 + 2 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
dev.off()

#### sensitivity analysis broken down by educational level ####
online.lowedu.data <- subset(online.data, college == 0)
online.highedu.data <- subset(online.data, college == 1)
placebo.lowedu.data <- subset(placebo.data, college == 0)
placebo.highedu.data <- subset(placebo.data, college == 1)

## high education
DQ.est.online.highedu <- mean(online.highedu.data$DQ, na.rm = TRUE)
DQ.se.online.highedu <- sqrt(DQ.est.online.highedu * (1 - DQ.est.online.highedu) / 
                            sum(! is.na(online.highedu.data$DQ)))

list.result.online.highedu <- lm_robust(list ~ treatment, data = online.highedu.data, 
                                        fixed_effects = ~ block)
list.est.online.highedu <- as.numeric(list.result.online.highedu$coefficients)
list.se.online.highedu <- as.numeric(list.result.online.highedu$std.error)

list.result.placebo.highedu <- lm_robust(list ~ treatment, data = placebo.highedu.data, 
                                         fixed_effects = ~ block)
list.est.placebo.highedu <- as.numeric(list.result.placebo.highedu$coefficients)
list.se.placebo.highedu <- as.numeric(list.result.placebo.highedu$std.error)

pi.M.highedu <- array(NA, c(76, 3, 6))
for (i in 1:3) {
  corrected.list.est.highedu <- as.numeric(list.est.online.highedu - 
                                             list.est.placebo.highedu) / (1 - lambda)
  corrected.list.se.highedu <- as.numeric(sqrt((list.se.online.highedu ^ 2 + 
                                                  list.se.placebo.highedu ^ 2) / 
                                                 (1 - lambda) ^ 2))
  corrected.SDB.se.highedu <- as.numeric(sqrt(DQ.se.online.highedu ^ 2 + 
                                                (list.se.online.highedu ^ 2 + 
                                                   list.se.placebo.highedu ^ 2) / 
                                                (1 - lambda) ^ 2))
  pi.M.highedu[, i, 1] <- pi.P[i] + corrected.list.est.highedu
  pi.M.highedu[, i, 2] <- pi.P[i] + corrected.list.est.highedu + 
    qnorm(0.025) * corrected.list.se.highedu
  pi.M.highedu[, i, 3] <- pi.P[i] + corrected.list.est.highedu + 
    qnorm(0.975) * corrected.list.se.highedu
  pi.M.highedu[, i, 4] <- DQ.est.online.highedu - 
    (pi.P[i] + corrected.list.est.highedu)
  pi.M.highedu[, i, 5] <- DQ.est.online.highedu - 
    (pi.P[i] + corrected.list.est.highedu) + 
    qnorm(0.025) * corrected.SDB.se.highedu
  pi.M.highedu[, i, 6] <- DQ.est.online.highedu - 
    (pi.P[i] + corrected.list.est.highedu) + 
    qnorm(0.975) * corrected.SDB.se.highedu
}

## low education
DQ.est.online.lowedu <- mean(online.lowedu.data$DQ, na.rm = TRUE)
DQ.se.online.lowedu <- sqrt(DQ.est.online.lowedu * (1 - DQ.est.online.lowedu) / 
                               sum(! is.na(online.lowedu.data$DQ)))

list.result.online.lowedu <- lm_robust(list ~ treatment, data = online.lowedu.data, 
                                        fixed_effects = ~ block)
list.est.online.lowedu <- as.numeric(list.result.online.lowedu$coefficients)
list.se.online.lowedu <- as.numeric(list.result.online.lowedu$std.error)

list.result.placebo.lowedu <- lm_robust(list ~ treatment, data = placebo.lowedu.data, 
                                         fixed_effects = ~ block)
list.est.placebo.lowedu <- as.numeric(list.result.placebo.lowedu$coefficients)
list.se.placebo.lowedu <- as.numeric(list.result.placebo.lowedu$std.error)

pi.M.lowedu <- array(NA, c(76, 3, 6))
for (i in 1:3) {
  corrected.list.est.lowedu <- as.numeric(list.est.online.lowedu - 
                                             list.est.placebo.lowedu) / (1 - lambda)
  corrected.list.se.lowedu <- as.numeric(sqrt((list.se.online.lowedu ^ 2 + 
                                                  list.se.placebo.lowedu ^ 2) / 
                                                 (1 - lambda) ^ 2))
  corrected.SDB.se.lowedu <- as.numeric(sqrt(DQ.se.online.lowedu ^ 2 + 
                                                (list.se.online.lowedu ^ 2 + 
                                                   list.se.placebo.lowedu ^ 2) / 
                                                (1 - lambda) ^ 2))
  pi.M.lowedu[, i, 1] <- pi.P[i] + corrected.list.est.lowedu
  pi.M.lowedu[, i, 2] <- pi.P[i] + corrected.list.est.lowedu + 
    qnorm(0.025) * corrected.list.se.lowedu
  pi.M.lowedu[, i, 3] <- pi.P[i] + corrected.list.est.lowedu + 
    qnorm(0.975) * corrected.list.se.lowedu
  pi.M.lowedu[, i, 4] <- DQ.est.online.lowedu - 
    (pi.P[i] + corrected.list.est.lowedu)
  pi.M.lowedu[, i, 5] <- DQ.est.online.lowedu - 
    (pi.P[i] + corrected.list.est.lowedu) + 
    qnorm(0.025) * corrected.SDB.se.lowedu
  pi.M.lowedu[, i, 6] <- DQ.est.online.lowedu - 
    (pi.P[i] + corrected.list.est.lowedu) + 
    qnorm(0.975) * corrected.SDB.se.lowedu
}

#### Figure A8 ####
cairo_pdf("Figure_A8a.pdf", width = 6, height = 4, pointsize = 8)
layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(4, 4, 1, 1), oma = c(0, 0, 2, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.74, 0.01), seq(0.74, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 2], 
          rev(pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.highedu[, 1, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.72, 0.01), seq(0.72, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 2], 
          rev(pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.highedu[, 2, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.71, 0.01), seq(0.71, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 2], 
          rev(pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.highedu[, 3, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.74, 0.01), seq(0.74, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 5], 
          rev(pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.highedu[, 1, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 1, 1] < 1, 1, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.72, 0.01), seq(0.72, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 5], 
          rev(pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.highedu[, 2, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 2, 1] < 1, 2, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.71, 0.01), seq(0.71, 0.24, -0.01)), 
        c(pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 5], 
          rev(pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.highedu[, 3, 1] < 1], 
      pi.M.highedu[pi.M.highedu[, 3, 1] < 1, 3, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Approval rate of \"Paul Nueva\" = 0", at = 0.19, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.05", at = 0.19 + 1 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.1", at = 0.19 + 2 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
dev.off()

cairo_pdf("Figure_A8b.pdf",
          width = 6, height = 4, pointsize = 8)
layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(4, 4, 1, 1), oma = c(0, 0, 2, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.65, 0.01), seq(0.65, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 2], 
          rev(pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.lowedu[, 1, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.63, 0.01), seq(0.63, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 2], 
          rev(pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.lowedu[, 2, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.61, 0.01), seq(0.61, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 2], 
          rev(pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 3])), 
        border = NA, col = "gray80")
box()
abline(h = DQ.est.online, lty = 3)
lines(lambda[pi.M.lowedu[, 3, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 1])
text(0.22, DQ.est.online + 0.03, "DQ estimate", pos = 4)
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Approval rate of Duterte", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.65, 0.01), seq(0.65, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 5], 
          rev(pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.lowedu[, 1, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 1, 1] < 1, 1, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.63, 0.01), seq(0.63, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 5], 
          rev(pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.lowedu[, 2, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 2, 1] < 1, 2, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.24, 0.71), ylim = c(-1, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.24, 0.61, 0.01), seq(0.61, 0.24, -0.01)), 
        c(pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 5], 
          rev(pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 6])), 
        border = NA, col = "gray80")
box()
abline(h = 0, lty = 3)
lines(lambda[pi.M.lowedu[, 3, 1] < 1], 
      pi.M.lowedu[pi.M.lowedu[, 3, 1] < 1, 3, 4])
mtext("Proportion of non-strategic misreporters", 
      side = 1, line = 2.5, cex = 0.7)
mtext("Social desiability bais", 
      side = 2, line = 2.5, cex = 0.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Approval rate of \"Paul Nueva\" = 0", at = 0.19, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.05", at = 0.19 + 1 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
mtext("Approval rate of \"Paul Nueva\" = 0.1", at = 0.19 + 2 / 3, line = 0.5, 
      cex = 0.8, outer = TRUE)
dev.off()